% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Replication "Deconstructing the Yield Curve"
% Crump and Gospodinov (2024)
% Date: 26-JUL-2024
% Function for bootstrapping VAR models
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
function bsOut = bootstrapVAR(Y1,Y2,B,B_BC_YLD,B_BC_EXT,doBH,biasCorrectExplosive)

    T = size(Y1,1);
    p1 = size(Y1,2);            
    varY1 = var1BC(Y1,B_BC_YLD,biasCorrectExplosive);
    tmpV1 = varY1.VBC;
    bsY1 = NaN(T,p1,B);
    % Need to demean the VAR residuals because using BC (no effect when B_BC==0)
    V1 = tmpV1 - ones(T,1)*mean(tmpV1, 'omitnan');
    %
    if ~isempty(Y2)
        p2 = size(Y2,2); 
        varY2 = var1BC(Y2,B_BC_EXT,biasCorrectExplosive);    
        tmpV2 = varY2.VBC;
        bsY2 = NaN(T,p2,B);
        % Need to demean the VAR residuals because using BC (no effect when B_BC==0)
        V2 = tmpV2 - ones(T,1)*mean(tmpV2, 'omitnan');
    end
    %
    for b = 1:B
        if ~doBH            
            V1b = mvnrnd(zeros(size(V1(2:end,:))),varY1.SigBC);            
            tmpBsY1 = buildVAR(varY1.muBC, varY1.PhiBC, V1b', Y1(1,:)');
            if ~isempty(Y2)
                V2b = mvnrnd(zeros(size(V2(2:end,:))),varY2.SigBC);
                tmpBsY2 = buildVAR(varY2.muBC, varY2.PhiBC, V2b', Y2(1,:)');
            end                        
        else
            idxRand = randi(T-1,[T-1,1])+1;
            tmpBsY1 = buildVAR(varY1.muBC, varY1.PhiBC, V1(idxRand,:)', Y1(1,:)');
            if ~isempty(Y2)
                tmpBsY2 = buildVAR(varY2.muBC, varY2.PhiBC, V2(idxRand,:)', Y2(1,:)');
            end
            %
        end
        bsY1(:,:,b) = [Y1(1,:)' tmpBsY1]';
        if ~isempty(Y2), bsY2(:,:,b) = [Y2(1,:)' tmpBsY2]'; end        
    end
    %
    bsOut.bsY1 = bsY1;
    if ~isempty(Y2), bsOut.bsY2 = bsY2; end
end